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Gaussian Quantum Monte Carlo (GQMC) is a stochastic phase space method for fermions with 
positive weights. In the example of the Hubbard model close to half filling it fails to reproduce all 
the symmetries of the ground state leading to systematic errors at low temperatures. In a previous 
work [Phys. Rev. B 72, 224518 (2005)] we proposed to restore the broken symmetries by projecting 
the density matrix obtained from the simulation onto the ground state symmetry sector. For ground 
state properties, the accuracy of this method depends on a large overlap between the GQMC and 
exact density matrices. Thus, the method is not rigorously exact. We present the limits of the 
approach by a systematic study of the method for 2 and 3 leg Hubbard ladders for different fillings 
and on-site repulsion strengths. We show several indications that the systematic errors stem from 
non- vanishing boundary terms in the partial integration step in the derivation of the Fokker-Planck 
equation. Checking for spiking trajectories and slow decaying probability distributions provides 
an important test of the reliability of the results. Possible solutions to avoid boundary terms 
are discussed. Furthermore we compare results obtained from two different sampling methods: 
Reconfiguration of walkers and the Metropolis algorithm. 

PACS numbers: 71.27.+a, 71.10.-w, 71.10.Fd 71.15.-m 



I. INTRODUCTION 

One of the biggest unresolved problems in computa- 
tional physics is the negative sign problem for fermionic 
and frustrated systems. Although it is not possible to 
solve it in general 0] there is hope to find solutions 
for specific models. Gaussian Quantum Monte Carlo 
(GQMC) |,i claims to be a sign-free ab-initio method 
for the general electronic structure problem. Pirst results 
for the Hubbard model looked very promising. However, 
in a previous paper [l| we have shown that in the vicinity 
of half filling systematic errors in the energy and other 
quantities occur, and that the method fails to reproduce 
the symmetries (SU(2) spin, translation, and lattice sym- 
metries) of the Hamiltonian. This broken symmetry can 
be restored by projecting the low temperature density 
matrix from the simulation onto the ground state sym- 
metry sector, such that excitations from other sectors 
are projected out. We have found that observables eval- 
uated with the projected density matrix agree with exact 
ground state results, as we have shown for Hubbard mod- 
els up to system sizes of 6x6. 

One important aim of this paper is to investigate the 
origin of the symmetry breaking and the systematic er- 
rors. Our results suggest that they appear due to non- 
vanishing boundary terms in the partial integration step 
in the derivation of the Pokker-Planck equation. We show 
that changing the sampling method does not help to 
avoid this problem. Due to the inherent inaccuracy of 
the density matrix as produced by the GQMC, it is not 



clear that symmetry projection schemes will produce cor- 
rect ground state properties. It is therefore important to 
analyze the limits of this method. To this end we present 
a systematic study of Hubbard ladders. 

The paper is organized as follows: The next section 
provides a summary of the method and the derivation of 
the stochastic differential equations (SDEs). A more de- 
tailed description was previously given in [5| . Section IIIII 
addresses the origin of the systematic errors. We show 
that slow decaying power law tails in probability distri- 
butions can cause two different kinds of problems. First, 
they may lead to diverging variances of observables, mak- 
ing a Monte Carlo sampling useless. Second, boundary 
terms may appear in the partial integration step in the 
derivation of the Pokker-Planck equation. In the pres- 
ence of boundary terms the SDEs are no longer valid, 
and by neglecting them a systematic error is introduced. 
This problem has been encountered before in the context 
of stochastic phase space methods for bosonic systems 
[1] and was solved for specific models with the help of 
stochastic gauges ■ A side effect of boundary terms is 
the presence of spiking trajectories, therefore checking for 
spikes in the phase space variables is an important test 
of the reliability of the results. In section HVl we discuss 
results from the Metropolis algorithm which leads to the 
same systematic errors as the reconfiguration scheme of 
walkers [1] we usually use. 

Empirically we have seen that one of the major conse- 
quences of fat tailed distributions shows up in the viola- 
tions of symmetries. Hence, imposing symmetry projec- 
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tions prior to measurements can potentially correct for 
sources of systematic errors. In section |V] we present a 
systematic study of the GQMC method with symmetry 
projection (PGQMC) for 2 and 3 leg Hubbard ladders 
for different fillings and on-site repulsion strengths, and 
compare the results with Density Matrix Renormaliza- 
tion Group (DMRG) calculations. Symmetry pro- 

jection is successful in removing systematic errors in all 
cases where the overlap of the density matrix with the 
ground state symmetry sector is not too small. However, 
for a small overlap systematic errors may still be present. 
In the outlook in section IVII we refer to recent promis- 
ing improvements of the projection method by Aimi and 
Imada fnl. 



II. SUMMARY OF THE METHOD 

Let us briefly recall the derivation of the SDEs which is 
a standard procedure for various stochastic phase space 
methods. A more detailed derivation can be found in Ref. 
0, 0] ■ The starting point is an expansion of the system 
density operator in an ovcr-complctc operator basis 



P{r) 



dXP{X,T)A{X), 



(1) 



where r is the inverse temperature and the probability 
density P is normalized / dAP(A, r) = 1. The A(A) are 
the Gaussian operator basis elements of the normal or- 
dered form 



A(n,0) = Sldct(l - n) 



r2+(n-- 



(2) 



with fct (respectively c) being an Ng dimensional vec- 
tor of creation (respectively annihilation) operators, n is 
an Ns X Ns real matrix of phase space variables and Ng 
denotes the number of states. det(l — n) is the normal- 
ization term such that Tr[A(n, il)] = il. Thus, ^l plays 
the role of a weighting factor. 

The imaginary time evolution of the density operator 
reads: 



dT 



P{r) 



(3) 



Introducing the expansion ((T|) for p leads to 



d 

d^ 



dXP{X,T)A{\) = 



dXP{X,T) H,A{X) 



(4) 

With the help of differential properties of the operator 
basis derived in Ref. Q the action of the Hamiltonian 
on the operator basis element can be transformed into 
an operator L containing first and second order deriva- 
tives with respect to the phase space variables A, and we 
formally write 



rfA— P(A,t)A(A) = / dAP(A,T)L[A(A)]. (5) 



In the next step we perform a partial integration where 
we take all derivatives acting on the basis element in front 
of P, which is denoted by the new operator L': 



dXL'[P{X,T)]A{X) + boimdary terms. 



(6) 



Depending on the nature of the distribution P, bound- 
ary terms from the partial integration arise. Let us first 
assume that these boundary terms vanish, so that we can 
compare integrands on both sides to obtain 



-P(A,r)=i'[P(A,T)], 



(7) 



where we have omitted the Gaussian basis element on 
both sides. This corresponds to a Fokker-Planck equation 
describing the evolution of the distribution function P in 
(imaginary) time. If L' is of the form 



L' 



af}k 



d_ 



d 



(8) 



with Aa and P„ real functions wc can derive real valued 
(Stratonovich) SDEs 



dK{T) = A^{X)dT + V Bi{X)dWk{T 



k 



(9) 



with Wiener increments dWk{T) defined by the correla- 
tions {dWkdWk') = drSkk' and the mean {dWkXr)) = 0. 
The explicit forms of the functions ^^.(A) and B^{X) 
for the Hubbard model can be found in appendix [X] 
The form of L' is not unique but can be modified by 
gauge degrees of freedom. In Ref. [1] the "Fermi gauge" 
nfj^ — fiii„ = is used to obtain real valued SDEs with 
positive weights. Adding such terms clearly does not 
modify the expectation value of the Hamiltonian, but 
changes the resulting Fokker-Planck equation. 



HI. SOURCES OF SYSTEMATIC ERRORS 

A. Systematic errors in the Hubbard model 

We have tested the GQMC method for the Hubbard 
model given by the Hamiltonian 



H 



{t.j),a 



(10) 



with nearest neighbor hopping strength t, on-site repul- 
sion U and chemical potential fi. The corresponding 
stochastic differential equations derived under the as- 
sumption of vanishing boundary terms can be found in 
appendix |^ 

As already pointed out in Ref. [l[ the GQMC method 
works well for small electron interaction U jt and away 
from half filling (Fig. [1] upper plot). In this regime the 
ground state is very well described by a paramagnetic 
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FIG. 1: (Color online) Energy as a function of inverse tem- 
perature P obtained from exact diagonalization (dashed lines) 
and GQMC runs (solid lines) for a 2 x 2 Hubbard model. In 
the upper plot the system is far away from half filling: (7 = 1, 
t = 1, fj, — —1, averaged over 40,000 trajectories. The exact 
result is reproduced within the statistical error bars. In the 
lower plot, close to half filling, the energy from the simulation 
is systematically too high {U = 4, t — 1, /j, = 1, averaged over 
480,000 trajectories). 



mean field solution which is exactly reproduced by the 
GQMC approach. Close to half filling and with a big 
on-site repulsion the simulation results exhibit system- 
atic errors (Fig. [1] lower plot). The energy agrees with 
the exact result down to a certain temperature, but for 
lower temperatures the mean energy is systematically too 
high. In Ref. [l| it was found that the solution of the 
simulation does not preserve SU(2) spin rotation symme- 
try. This gave the motivation to develop the projection 
scheme as described in section |Vl But the reason for this 
symmetry breaking has not been found so far. In sec- 
tion [TTrC] we suggest that the systematic errors and the 
symmetry breaking stem from non-vanishing boundary 
terms. 



B. Power law tails in the probability distribution 

of observables 



In this section we discuss how a power law tail in the 
probability distribution of an observable X can lead to 
problems in the Monte Carlo sampling. The error on 
the expectation value {X) obtained from a Monte Carlo 
simulation is given by AX/ \fM where M is the number 
of independent samples and AX the variance 



AX = {{X^) {Xf) 



1/2 



(11) 



If the variance of X diverges, then also the error bar 
of our Monte Carlo result diverges. Thus, to obtain a 
meaningful result from the sampling, the variance has 
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FIG. 2: (Color online) Log-log plot of the distribution of the 
transverse spin susceptibility for the 2x2 Hubbard model 
at half filling with U/t = 4 at low temperatures (/? = 10). 
The slope of the power law tail yields p « 2.6, therefore the 
variance is not defined. 



to be well defined. The m**^ moment of the probability 
distribution P{X) is given by 



{X' 



j X"'P{X)dX. 



(12) 



If the probability distribution exhibits a power law tail 
P{X) cx X~P for X oo then only moments m < p — 1 
of P{X) converge, because 

J X"'P{X)dX ^ J X"'X-PdX = 

X("^-p)dX oo, form-p>-l. (13) 



Therefore, to obtain a finite mean corresponding to the 
first moment (m = 1) of P{X), p has to be bigger than 2. 
For a finite variance the integral has to converge also for 
TO = 2, which inquires p > 3. If p < 3 we do not obtain 
a meaningful result from a Monte Carlo sampling. 

We have found a diverging variance for the transverse 
spin susceptibility Xxy a-t low temperatures in a simula- 
tion exhibiting systematic errors. The slope of the distri- 
bution P{xxy) in the log-log plot in Fig. [^lyields p sa 2.6 
such that the second moment is not defined. We also 
observe systematic errors in the energy. However, the 
variance of the energy is always well defined in our simu- 
lations. Thus, the systematic errors found in the energy 
cannot be explained by an ill defined Monte Carlo sam- 
pling. In the next section we demonstrate that power law 
tails in the distribution of the phase space variables can 
cause a different kind of problem, namely non-vanishing 
boundary terms. 



C. Non- vanishing boundary terms 

Boundary terms (BTs) from the partial integration in 
Eq. [5] appear if certain moments of the high-dimensional 
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FIG. 3: Appearance of a power law tail in the distribution 
P(r) (log- log plot) at low temperatures (large The sim- 
ulation parameters are U — 4, t = 1, fi — 1 and 480,000 
trajectories. 



are the implications of BTs on our simulation? We have 
derived SDEs under the assumption of non-vanishing 
BTs. As soon as they become non- negligible, the SDEs 
are strictly speaking no longer valid, such that we obtain 
a wrong distribution P{X,t) for t > Pbt- We suggest 
that this is the origin of the systematic errors in the en- 
ergy and in other quantities. The problem of boundary 
terms has been discussed in detail for several bosonic sys- 
tems 01 ■ The following analysis is done in a similar 
spirit as for these systems. 

A first test for the presence of BTs is to measure 
the radial averaged distribution P{r) ^ r^^ with r = 

\j^ijG '^ijtr different inverse temperatures [3. The 
power law tail in Fig. [3] reaches a slope of p « 3.8. This 
already indicates that power law tails are also present in 
the high-dimensional distribution P(A) of all phase space 
variables. From the partial integration step in Eq. [Hlwe 
can find the explicit expressions for the BTs. They are 
of the form 



probability distribution P(A, r) do not converge. Ini- 
tially, at infinite temperature, P{\, r = 0) is a delta- 
function, for which all moments converge and therefore 
there are no BTs at the beginning. As the distribution 
P evolves in imaginary time it becomes broader and we 
show below, that it develops slow decaying power law 
tails, such that BTs cannot be excluded anymore. Typi- 
cally BTs appear after a specific imaginary time r > (3bt- 
But other situations are known 0], where the BTs are 
present only for a short time and disappear again. What 



)M(.)(A("V(A)A(A) 



(14) 



Boundary 



where a enumerates the different BTs, M^^) depends on 
the phase space variables up to fourth order and the inte- 
gral is taken over all phase space variables except the one 
the partial integration has been carried out for. The basis 
element A(A) also depends on the phase space variables. 
It includes terms in riyo- up to 2Ns th order. We present 
here one specific example of a non-vanishing boundary 
term for a 2 site model at half filling with U /t = 100: 



lim lim ... lim / rffiiij / dn 



ailT 



ai3T 



ai4T 

13T / (irii4|. 

ai4T 



driNNi n52T^(A) 



+oi2-r 



-ai2t 



(15) 



r 



Comparing with Eq. [TJ] this term corresponds to 
M(q,)(A''"') = ni2\ni2-\, and from the expansion of A(A) 
we took the term proportional to ni2-\ , leading to an in- 
tegrand which is cubic in the phase space variable ni2|. 
Eq. [15] stems from a partial integration with respect to 
the variable ni2| . The remaining integral is taken over all 
variables nyg- ^ "■12? i yielding the marginal distribution 



P{ni2'\)- 



lim "?2T-P("i2T)|tI''! 



The fit to the power law tail of the distribution P(ni2|) 
in Fig. m yields an exponent of p k, 2.6. Therefore, with 
~ '^12'!^ this term does not vanish and we there- 
fore cannot neglect it in the partial integration step. This 
BT is particularly simple to analyze because it involves 
only a one-dimensional distribution. For other BTs with 



mixed variables one would need to study a distribution of 
several variables. By studying all possible terms appear- 
ing one could find the minimal exponent pbt necessary 
to exclude BTs. But as discussed in Ref. there are 
simpler ways to detect BTs as we show in the next sec- 
tion. 

The question may arise why power laws actually occur. 
The problem lies in the strong multiplicative noise term 
in our SDEs. The noise is amplified by the phase space 

(16) 

variables riyo- themselves (the diffusion term in Eq. [9] 
is even quadratic in nij„). This has strong consequences 
on the functional form of the distribution P{X). In Ref. 

it is shown that the smallest multiplicative noise in 
the linear Langevin equation changes the Gaussian sta- 
tionary distribution to one with a power-law tail. Power 
laws can arise from multiplicative stochastic processes as 
naturally as Gaussian distributions from processes with 
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FIG. 4: (Color online) Power law tail in the distribution 
^'('^121) = P{ni2i) for 2 site system at half filling for U/t = 
100, /3 = 3, 800,000 walkers. 



additive noise [iSj . The question remains if at finite time 
r the distribution will always have a finite cut-off at a 
certain distance such that the BTs would always van- 
ish. However, we have not found such a cut-off in our 
simulations. 



D. Spiking trajectories 

There are other indications that the boundary terms in 
equation ([6|) do not vanish. According to Ref. [a| spiking 
behavior of the trajectories imply that BTs are likely to 
be significant. Such near-singular trajectories do large 
excursions in phase space for a very short time (within a 
few time steps). Spikes could also stem from an unstable 
integration scheme. It is therefore important to use a 
stable integrator. Fig. [5] shows an example of a sharp 
spike in the energy. Such extreme trajectories lead also 
to a sudden increase of the statistical error of observables. 

As mentioned in the last section, there are no BTs at 
the beginning of the simulation, as we start from a delta 
function for the distribution P{\,t = 0). They appear 
at a specific inverse temperature Pbt as we integrate the 
SDEs towards lower temperatures. The first appearance 
of a spike should give an estimate of (3bt- In the example 
of Fig. [1] (lower plot) the first spiking walker shows up for 
f3BT ~ 1-5, which is in good agreement with the inverse 
temperature, at which the energy starts to deviate from 
the exact result. 

For small interaction U /t, where we obtain correct 
results, not a single spike can be observed. Therefore 
testing for spikes provides a good indicator, whether the 
GQMC results are reliable. 
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FIG. 5: (Color online) Upper figure: Mean energy from the 
simulation compared with the exact energy. The spike at 
/3 ~ 5.24 is caused by a single extreme trajectory shown in 
the lower figure (solid line). The spike occurs within a few 
time steps. The simulation is done at half filling for a 2 x 2 
lattice with U/t = 4, and 10,000 trajectories. 



are possible. Thanks to the overcompletness of the basis 
several solutions of the distribution P(A, t) exist. With 
appropriate gauges the boundary terms could possibly 
be removed , leading to a more compact distribution. 
The so-called drift gauges could be used to avoid nearly- 
diverging trajectories which cause power law tails. The 
tradeoff is to introduce noise into the equation for the 
weight. Stochastic gauges have been successfully applied 
for several models [7| . A future analysis will show if sim- 
ilar techniques can be applied to the Hubbard model to 
solve the present problems. 



IV. METROPOLIS ALGORITHM 



The equation for the weight Q{t) (equation (jA2p from 
the appendix) of a walker can be integrated. The weight 
as a function of time then reads 



n{P) = exp 



H{n{T)) dr 



(17) 



E. Stochastic gauges 

As already mentioned for a given Hamiltonian the 
SDEs are not unique, but different choices of "gauges" 



The weight and the variance of the weight thus grow ex- 
ponentially, yielding the need of some importance sam- 
pling. We usually use a reconfiguration scheme similar to 
the one used in the Green function Monte Carlo method 
Q . A further way to sample the distribution is to use the 
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Metropolis - Hasting algorithm which was recently pro- 
prosed by Dowling et al. • In this section we want to 
briefly summarize the basic ideas of the algorithm within 
the framework of the GQMC and then present some re- 
sults. 



consider A'' time steps. Then one sample is given by a 
noise vector w G R*^^ and the phase space variables for 
the final time t = N Ar can be considered as a function 
of the noise vector, i.e. n(iS) and ^{w). w is a normal 
distributed random variable with a probability density 
given by 



A. Metropolis Algorithm 

Starting from an arbitary state s„ one uses a suitable 
candidate generation function q{s^ s') in order to create 
a proposed step. This step is then accepted s„+i = s' 
with the probability 



' TT{s')q{s\s) 

a = mm | , 1 

7r(sjg(s, s ) 



(18) 



If the move is rejected, the old state is kept Sn+i = Sn- 
One can prove that this algorithm fulfills the detailed 
balance so the resulting chain samples the target density 
TT correctly [l^. 
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P{w) 



cxp 



w 



(19) 



The expectation value of an observable O can thus be 
written as 



O) ^ 



Trp 6 _ {n{w) 0{w)) 

Trp ^ {^{w)) 
/ P{w) Q{w) Oiw) d'^'^i 
/ P{w) n{w) d^^^w 



(20) 



Now, one can apply the Metropolis-Algorithm to create 
samples (wi); which are distributed like 



tt{v3) = P{w) n{v3), 



(21) 



and finally the expectation value of the observable is 



given by < O >= -^^i^iO{'Wi). It is convenient to 
use a candidate generating density which obeys 

(l{w,w') ^ P{w') 

q{w\ w) P{w) ■ ^ ' 

The acceptance rate is then just the quotient of the 
weights ri, i.e. 



a 



(23) 



We used a simple candidate generation function which al- 
ters each component of the noise vector with a fixed prob- 
ability r, i.e. drawing approximately rMN new Gaus- 
sian numbers. One can easily adapt this algorithm to an 
adaptive step size (see appendix [B]) using a dynamical 
enlarged noise space. 
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FIG. 6: (Color online) Probability distributions of the total 

^ 2 

spin squared S (a) and the energy E (b) evaluated at /3 = 20 
obtained by using the Metropolis algorithm and the reconfig- 
uration scheme (2x2 Hubbard model, U/t — 4, half filled, 
100,000 samples each) 

The usual way to solve SDEs (as in Eq. ([9])) is to apply 
an Euler - Marujama method (either implicit or explicit). 
For each time step one needs a fixed number of Gaussian 
distributed random numbers. Let this number be M and 



TABLE I: Comparison between the Metropolis and the re- 
configuration results, (2x2 Hubbard model at half filling, 
U/t = 4, evaluated at /3 = 20, 100,000 samples). The error 
bar on 5*^ is ill defined because the variance of 5*^ diverges. 



B. Results 

Now we want to present some results using the 
Metropolis algorithm. The system which we discuss here 
is a 2 X 2 Hubbard model at half filling with U/t = 4. The 
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behaviour of this model is representative of that seen in 
other systems, and seems to be quite generic. The SDEs 
are solved using an implicit Euler scheme with an adap- 
tive time step with Armax = 5 • 10~* (see appendix [B]) 

The data using the reconfiguration scheme is obtained 
from 100,000 walkers and by applying the scheme ev- 
ery Ar^c = 0.05. For the Metropolis algorithm a typical 
chain length is 1,000 after a burn- in time of the order of 
100 Metropolis steps. The Metropolis algorithm has one 
major drawback compared to the reconfiguration scheme, 
namely one needs to fix a specific target time (}. When 
using the reconfiguration, one is in principle able to ob- 
tain values for any intermediate time and not just for 
the final one. Therefore the computational effort to cal- 
culate the observables for all times j3 = [0, (3finai\ is 
much bigger for the Metropolis algorithm than for the 
reconfiguration. An interesting variant of the Metropolis 
algorithm which allows the calculation of observables for 
intermediate times is presented in Ref. [ll|. 

Table [J shows the expectation values of the total spin 
squared and H at the target time /3 = 20. The ob- 
servable is chosen since after some time the variance 
of this observable is not finite anymore, thus yielding a 
good test whether changing the sampling method im- 
proves the results. For the results for the Metropolis 
algorithm 100 Markov chains with each 1,000 steps after 
the burn-in time are used, thus yielding 100,000 samples, 
which is the same number as the one used for the recon- 
figuration scheme. One clearly sees that the results do 
not change significantly, the energy is slightly improved 
but S'^ gets worse. To further investigate the effect of the 
sampling method, the probability distributions of the two 
observables were also calculated (see Fig. [5]). Again both 
methods produces almost the same probability distribu- 
tion, especially the slow decaying power law tails of S"^ 
are still present. Using the Metropolis sampling therefore 
does not seem to change the results significantly. 



V. SIMULATION OF HUBBARD LADDERS 

In Ref. [H we have found that GQMC fails to repro- 
duce all the symmetries of the Hamiltonian at low tem- 
peratures. We proposed to restore the broken symme- 
tries by projecting the density matrix obtained from the 
simulation psim onto the ground state symmetry sector, 

Pproj = PpsimP\ (24) 

where P is the corresponding projection operator. For 
example projecting the density matrix onto the S = 
sector filters out spin excitations and restores the SU(2) 
rotation symmetry in spin space. 

The discussion in section Hill suggests that the symme- 
try breaking is directly related to the presence of bound- 
ary terms. In Ref. [l| we obtained correct results for 
the ground state by an appropriate symmetry projec- 



tion, which implies that even in the presence of bound- 
ary terms, the correct ground state is still included in the 
density matrix, but mixed with excited states, which we 
can project out. The question remains, if such a projec- 
tion can always be done. The aim of the current section 
is to show the limits of this projection method (PGQMC) 
for the example of Hubbard ladders for various lengths, 
interaction strengths and doping. Energy and correla- 
tion functions are compared with calculations from the 
Density Matrix Renormalization Group (DMRG [ol. [lOjl 
method which provides high precision results for quasi 
one dimensional systems. The DMRG results are cal- 
culated in a matrix product state basis using both the 
SU(2) spin and the SU(2) pseudospin symmetry flT^ . 

The PGQMC simulations in this section are done 
with 8,000 - 32,000 walkers, an adaptive time step with 
Arniax = 5 • 10^** (see appendix [B|) and an explicit Euler 
integration scheme. For some examples crosschecks have 
been made with more walkers, with an implicit Euler 
scheme and different quantization axis. The error bars 
stem from averaging over several projections at differ- 
ent imaginary times. They do not take into account the 
discretization error in the symmetry projection. For the 
comparison with the exact values we check if they are 
within two standard deviations {2a). If they lie outside 
we have to assume that besides the statistical error, there 
is also a systematic error present. 

Projection onto the ground state is only possible if 
there is a finite overlap (Tr[/5pi-oj]/Tr[psini]) between the 
density matrix and the ground state symmetry sector, or 
in other words, if the density matrix from the simulation 
contains a finite contribution of the ground state, which 
we can extract by the projection. The projection is al- 
ways made onto the S = sector and onto all possible 
momentum and parity sectors. Note that we only have 
translational symmetry along the x-axis, whereas along 
the y-axis we can distinguish odd and even parity. The 
ground state sector is the one with lowest energy, and 
we found that this sector always has the biggest overlap 
with the density matrix. 



A. Two leg ladders 

1. L—4 and varying U 

As already pointed out GQMC works well in the weak 
interacting case (small U/t) and systematic errors ap- 
pear for large U/t. For U/t = 1 the GQMC simu- 
lation results for the spin-spin and charge-charge cor- 
relation functions in Fig. [7] agree with the DMRG re- 
sults even without symmetry projection. Also the en- 
ergy EcQMC = —10. 117 ± 0.001 is correct compared to 
Edmrg = —10.118. 

For U/t = 2 we also obtain correct results without 
projection. Systematic deviations of the order 5% appear 
for U/t = 4, which are corrected by symmetry projection. 
As expected we observe an increase of the systematic 
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FIG. 7: (Color online) Spin-spin and charge-charge correla- 
tion function of the half filled 4x2 Hubbard ladder for U/t=l 
(without symmetry projection) showing perfect agreement 
with DMRG. The deviation from the DMRG result (dashed 
line) is shown in the inset. 



FIG. 9: (Color online) Energy of the half filled 4x2 Hubbard 
ladder in dependence of U/t. The deviation from the DMRG 
result (dashed line) is shown in the inset. Symmetry projec- 
tion corrects the systematic deviations in the energy from the 
GQMC simulation. 



deviations with increasing U/t. For U/t — 16 symmetry 
projection yields the correct energy but fails to reproduce 
all the spin-spin correlations at large distance correctly 
(see Fig. ID). In this case the average overlap of 28% 
is rather small. Thus symmetry projection yields better 
results for intermediate U/t but for large U/t > 8 a. small 
systematic error is still present. In Fig.[5]we have plotted 
the dependence of the energy on U/t. Without projection 
the systematic error grows with increasing U/t, whereas 
the results after projection agree with the exact result for 
all U/t. 



2. U/t — 4 and varying L 

111 this set of simulations we fixed the interaction 
strength to an intermediate value U/t = 4 and varied 
the system length L. For L > 4 we observe that the 
energy from the GQMC simulation is systematically too 
high (Fig. (TU]). The deviations are of order 2%. After 
symmetry projection the results are within 2a for sys- 
tem sizes up to L = 16. 

Excellent results for the correlation functions arc ob- 
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FIG. 8: (Color online) Spin-spin and charge-charge corre- 
lation function of the half filled 4x2 Hubbard ladder for 
U/t=16 after projection. The deviation from the DMRG re- 
sult (dashed line) is shown in the inset. A small systematic 
deviation is still present in the spin-spin correlations at large 
distance. 



FIG. 10: (Color online) Energy density of the half filled two 
leg Hubbard ladder with U/t — 4 in dependence of L. The 
deviation from the DMRG result (dashed line) is shown in 
the inset. The energy from GQMC is systematically too high 
(dots). The results are correct after symmetry projection for 
L < 16. 
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FIG. 11: Spin-spin correlation function of the half filled 8x2 
Hubbard model with U /t = 4. The PGQMC (dots) results 
agree with the DMRG result (dashed line). The deviation 
from the DMRG result is shown in the inset. 



tained with the PGQMC method for L = 8 (Fig. [TTl) 
and L = 10. Without projection the results are qual- 
itatively good, systematic deviations are of order 10% 
for the spin-spin correlations and of order 0.2% for the 
charge-charge correlations. For L = 16 (Fig. [T^ the val- 
ues of the spin-spin correlations at large distances tend to 
be too large (in absolute value), such that systematic er- 
rors may still be present. The problem is that the overlap 
with the ground state sector decreases with increasing L. 
For L = 4 it is typically of order 70% whereas for i = 10 
we find an overlap around 20%. For L > 16 it is only 
a few percent, such that the results from the projection 
method are not reliable anymore and we find systematic 
deviations even after projecting. 

At half filling the total number of particles stays con- 
stant, but the number of particles with spin up n-f is not 
necessarily equal to the number of particles with spin 
down n|. For big system sizes the simulation can get 
stuck in a configuration where 7^ ti^ . In this case 5"*°*"' 
is not equal to zero and therefore the overlap with the sec- 
tor 5"*°*"' = almost vanishes. A solution to this problem 
is to use the quantization axis along the x-direction, lead- 
ing to identical SDEs for and nj^ and therefore — n| 
is always guaranteed. However, even with this variant the 
results are not satisfying. The GQMC result shows big 
systematic deviations in the correlation functions at large 
distance. The overlap even becomes negative (and small 
in absolute value) for some projections, because many of 
the projected weights arc negative. This leads to cancel- 
lation between positive and negative weights, reminiscent 
of the sign problem in conventional QMC. 



3. Doped examples 

Next we present results for doped Hubbard ladders. 
The chemical potential ji is adjusted to obtain the desired 
number of electrons and we fix U/t = 4. The results 
for the 8x2 system with ritot = 14 are compatible with 
the DMRG results (Fig. [75]) . However, the error bars 
are bigger compared to the examples at half filling. The 
overlap with the ground state sector is only « 12%. 

We also find agreement for the 10x2 system with 
ntot = 18, even if the average overlap is only 6%. The en- 
ergy from DMRG -16.6393 lies within 2a of the PGQMC 
result -16.2244 ± 0.3206. Shghtly doped Hubbard mod- 
els are known to exhibit a strong sign problem in con- 
ventional QMC, therefore it is a considerable success to 
obtain correct results for this case. 



B. Three leg ladders 

Three leg Hubbard ladders are critical, thus the 
energy gap vanishes in the thermodynamic limit. As the 
low lying excitations lie closer to the ground state we 
expect that the density matrix from the simulation will 
contain more admixtures of excited states than for the 
two leg ladders. This would result in a smaller overlap 
with the ground state sector and thus a less efficient 
symmetry projection. We tested system sizes from 
L = 4toL=16. We obtain correct results up to L = 12 
(Fig. [H]) . For L = 16 the overlap becomes too small, 
leading to very large error bars (see Table |TT|. Note 
that in this case the energy obtained from GQMC alone 
is actually better than the energy afer the projection. 
However, the GQMC spin-spin correlations show large 
systematic deviations for large distances. Thus, neither 
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FIG. 12: Spin-spin correlation function of the half filled 16x2 
Hubbard model with U/t = 4. The PGQMC (dots) resuhs 
agree with the DMRG result (dashed line) for almost all dis- 
tances. The deviation from the DMRG result is shown in the 
inset. 
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of the two methods yield useful results for L = 16. 



L 


GQMC 


PGQMC 


DMRG 


Overlap 


4 


-9.042 ± 0.007 


-9.213 ± 0.015 


-9.2053 


57% 


6 


-13.54 ±0.05 


-13.69 ±0.10 


-13.7901 


54% 


8 


-17.98 ±0.02 


-17.98 ±0.23 


-18.2005 


43% 


12 


-26.54 ± 0.04 


-27.24 ±0.08 


-27.2717 


21% 


16 


-36.57 ±0.21 


-32.02 ± 2.93 


-36.2408 


2% 



TABLE II: Energies for the Lx3 Hubbard ladder at half filling 
The resuhs from PGQMC are within 2a up to i = 12. 



C. Summary of the Hubbard ladder results 

Let us summarize the results for the Hubbard ladders: 

• For small system sizes GQMC yields correct results 
for weak interaction {U/t < 2). Systematic devia- 
tions for intermediate interaction strength can be 
fixed by symmetry projection, but only for U jt not 
too large ^ jt < 8). For strong interaction we find 
systematic errors even for the symmetry projected 
result. 

• The overlap has to be big enough in order to get 
meaningful results. In our simulations overlaps of 
> 30% lead to correct results. Overlaps below 10% 
are definitely too small. We obtained some nice 
results with overlaps in between 10 — 30%, but the 
reliability is not guaranteed. 

• The overlap between the GQMC density matrix 
and the symmetry sector of the ground state de- 
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FIG. 13: Spin-spin correlation function of the doped 8x2 Hub- 
bard model with U/t = 4 and ritot = 14. The PGQMC (dots) 
results agree with the DMRG result (dashed line). The devi- 
ation from the DMRG result is shown in the inset. 
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FIG. 14: Spin-spin correlation function for the half filled 8x3 
Hubbard model with U/t — 4 after projection. The distances 
indicated by square brackets refer to distances on the middle 
leg of the ladder. The deviation from the DMRG result is 
shown in the inset. 



creases with increasing system size. In our exam- 
ples the overlap becomes too small for L > 16. Us- 
ing more walkers may help to increase the upper 
limit of L for the which PGQMC produces good 
values. 



VI. SUMMARY AND OUTLOOK 

The discussion in section UlIl suggests that the system- 
atic errors found in the Hubbard model close to half fill- 
ing stem from non-vanishing boundary terms from the 
partial integration step in the derivation of the SDKs. 
This problem has also been reported for bosonic systems 
0, and we observe similar side effects of the boundary 
terms, such as spiking trajectories. Thanks to the over- 
completeness of the Gaussian operator basis it is possi- 
ble to modify the SDEs without changing the physical 
system. The hope is to find appropriate gauges to ob- 
tain a faster decaying distribution function, such that 
the boundary terms vanish. The study of the Metropolis 
algorithm showed that it leads to the same systematic 
errors as the reconfiguration scheme of walkers. 

It is important to point out that for a large parameter 
range of the Hubbard model GQMC yields the correct 
results. Therefore it would be worthwhile applying the 
method also to other models. Checking for spikes and 
slow decaying probability distributions provides an im- 
portant test of the reliability of the results. 

We observed that the main effect of the boundary 
terms is that the solution from the simulation does not 
exhibit all the symmetries of the Hamiltonian. By pro- 
jecting the density matrix onto the ground state sym- 
metry sector it is possible to extract the correct ground 
state. In order to find the limits of the symmetry projec- 
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tion method wc have systematically tested it for Hubbard 
ladders and compared the results with DMRG calcula- 
tions. The overlap between the density matrix and the 
ground state sector has to be big enough (> 30%) in order 
to obtain good results. The results agree well for systems 
up to 32 sites and an on site repulsion U < 10, beyond 
these values the overlap becomes too small. However, we 
were able to obtain the correct values for doped ladders, 
which suffer from the negative sign problem in auxiliary 
field QMC. 

To conclude, even though GQMC is sign-free there are 
still unresolved problems. A future study will show how 
the boundary terms depend on the choice of gauges and 
if it is possible to avoid systematic errors even without 
projection. 

Recently Aimi and Imada (llj presented a new vari- 
ant of the projection method, called the pre-projection 
method, which allows to treat bigger system sizes. In- 
stead of projecting the density matrix after the simula- 
tion they incorporate the projection into the sampling, 
which leads to a better convergence towards the ground 
state. The price however is the occurence of negative 
weights, or in other words a sign problem which in many 
cases seems to be tractable. The results for Hubbard 
models up to system sizes 10x10 look very promising. For 
intermediate U the projected distributions decay much 
faster, so that no boundary terms seem to be present. As 
the method enables to simulate doped and frustrated sys- 
tems, it is one of the most promising ground state meth- 
ods for fermions currently available. However, further 
tests are needed to check the reliability of the method. 
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The Stratonovich stochastic differential equations for the 
Hubbard model read 

dn{T) = -QhdT (A2) 
dn^yp - A^,pdT + J2{KvpdW, + C'^,pdWl)iA3) 



with 



A,. 



5 EC 



X ^t(5(j,j) -I- \U\{nup - snu-p - ^)(5y + 



Buvp \l ^{nu^^n^^,^Sp^ ~ Snu^ln^vi5pi) 




CLp = \l ^{nui^rii^^Sp^ - snuiin,^\,5pi), (A4) 

where (i, j) denotes nearest neighbor pairs and the noise 
terms dWi are defined by the correlations {dWidWii) = 
drSiii and the mean {dWi{T)) = 0. We use the notation 

The drift term in the Ito formulation reads 



\Ito 

uvp 



'^ujp'^ivp ~^ '^ujp'^ivp) 



^ {^^{i-.j) ~ Unu-pSij + ^iSij) . (A5) 



APPENDIX B: ADAPTIVE TIME STEP 

To reduce the error from the time discretization of the 
SDEs we use an adaptive time step. Initially we choose a 
maximal step size Armax- Whenever any element of the 
drift term exceeds a certain threshold. 



max(A„„p • At) > u„ 

uvp 



(Bl) 



APPENDIX A: THE SDES FOR THE HUBBARD 
MODEL USED IN THE SIMULATION 

For the derivation of the real valued, positive weighted 



SDEs the "fermi gauge" [IJ nf, 
rewrite the interaction term as 



M 

2 



2 



was used to 



sign(C/). 



(Al) 



we divide the current time interval At by 2 and perform 
2 update steps with a reduced step size At/2. In each 
of these two steps the above condition is tested again 
and the step size is further divided by 2 if necessary. By 
proceeding in the same way for each reduced interval the 
step size can become arbitrarily small (limited by the 
machine precision). For example, the interval At may 
be split into 4 smaller intervals with step sizes At/2 -|- 
At/8 + At/8 + At/4. We usually choose u„iax between 
0.01 and 0.05. 
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